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presents a grand challenge and a major issue is the effective control of coherent 
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the task of determining the behaviour of the system. Here, we show how to de- 
termine open system Markovian dynamics of a quantum system with restricted 
initialisation and partial output state information. 
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1. Introduction 

Recently, much effort has been put into the design and realization of large scale 
quantum devices operating in the coherent regime. This has been spurred by the 
possibilities offered by quantum communication and information processing, from 
secure transmission, simulation of quantum dynamics, and the solution of currently 
intractable mathematical problems pQ. Many different physical systems have been 
proposed as basic architectures upon which to construct quantum devices, ranging 
from atoms, ions, photons, quantum dots and superconductors. For large scale 
commercial applications, it is likely that this will involve scalable engineered and 
constructed devices with tailored dynamics requiring precision control. 

Due to inevitable manufacturing tolerances and variation, each device will display 
different behaviour even though they may be nominally "identical". For their 
operation, they will need to be characterised as to their basic properties, response 
to control fields, and noise or decoherence [2]. We may also need to know how ideal 
is the system in the first place, for instance the effective Hilbert space; what we may 
assume to be a qubit may have dynamics involving more than two effective levels [3]. 
Extracting this information efficiently and robustly is crucial. 

In a laboratory setting, an experimentalist may have access to many tools with 
which to study a system, e.g. spectroscopy and external probes. In a production 
setting, provision of these extra resources may be difficult, expensive, or impossible 
to integrate with the device. It is therefore an important to understand what sort of 
characterisation can be performed simply using what is available in situ. Ideally, we 
would also like to be able to characterise the performance of a device with as little 
prior parameter data, e.g. how it responds to control fields, if this is the information 
which we are trying to obtain in the first place. Characterisation using only the in 
situ resources of state preparation and measurement, even where it is possible, is 
challenging, due to the increasing complexity of the signals, number of parameters 
to be estimated and the complexity of reconstructing a valid Hamiltonian from the 
resulting signal parameters. Robust and efficient methods of data gathering and 
analysis, preferably in as an automated way as possible, are therefore essential. 

System identification attempts to extract dynamics from minimal pre-existing 
system control resources. There is a trade-off between the universality of the dynamics 
that can be identified versus the resources available in terms of initialisation, access, 
and measurement. Here, we explore the issue of what can be extracted in practice, 
leaving aside the question of efficiency for future work. We show how we can extract 
the dynamics of an open quantum system with restricted state initialisation and 
only partial information of the evolving state. This relies on prior knowledge of the 
structure of the decoherence processes. Our results show that in some cases full state 
information without prior knowledge can lead to worse estimation than partial state 
information combined with limited a-priori knowledge about the expected model type. 

1.1. Quantum Process Tomography and System Identification 

In quantum process tomography (QPT) [3], what is determined is the discrete map 
between initial and final states. This requires in general the preparation of a complete 
set of input density operators, and quantum state tomography (QST) [5] of each of 
the output states, the result of which gives a completely positive map. However, 
this does not fully characterise the dynamical behaviour of the system such as the 



Open quantum system identification 



3 



Hamiltonian or Lindblad operators for Markovian systems. In principle, by conducting 
QPT at various times, it is possible to construct a time dependent set of Kraus 
operators [BJ, hence capturing the most general (not necessarily Markovian) open 
systems evolution but this is costly in both terms of initial state preparation and 
QST. We therefore look for methods which require fewer resources to determine the 
dynamics of a system, though perhaps at the expense of the ability to handle the most 
general evolutions, considering that for many systems the physics is sufficiently known 
so that the structure of the dynamics can be constrained, even if the precise values 
(which we are trying to determine) may not be known. 

Previous work has examined several system identification scenarios. Hamiltonian 
identification for a single qubit was introduced in [H [9] , Hamiltonian identification 
in the presence of decoherence was examined in [10j . extension to higher dimensional 
systems and the use of Bayesian/Maximum Likelihood techniques in [IT] and [12] . and 
model-based system parameter reconstruction in [T3]. Other work has also studied 
identification of system Hamiltonians for spin networks with restricted access [151 116j 
and limited system initialisation |17) . 

In the setting where preparation and measurement are limited to a single fixed 
basis it was found that signal and parameter complexity proved to be a major 
challenge even for relatively small systems and Hamiltonian evolution |11) . Including 
open system behaviour magnifies the problem even further. Besides increasing the 
number of parameters to be estimated, damping reduces the signal length and 
limits the amount of information that can be extracted [10]. To make the problem 
tractable we can try to incorporate prior information about the structure of the 
dynamics [JJ5] [T5]. In this case we find that even when only partial projective 
measurement data is available the signal may contain sufficient information to enable 
Hamiltonian reconstruction |13j . 

Here we examine problem of determining open quantum system dynamics based 
on restricted initial state preparation and incomplete state reconstruction, e.g., limited 
to measurements of a single observable. Using a combination of Bayesian/Maximum 
likelihood estimation of signal parameters [U] coupled with analysis of the equations 
of motion in the Laplace domain we study the problem of idcntifiability of models, 
comparing scenarios involving different types and amounts of information available for 
reconstruction. The results are applied to the problem of identifying the dynamics of 
a qubit in a Markovian environment to reconstruct both Hamiltonian and decoherent 
processes in various settings. 

2. Master equation for Markovian Case 

The evolution of an iV-level open quantum system in a Markovian environment is 
given by the dissipative master equation 

, . N 2 -l 



where p(t) is the state of the system at time t and if is a Hermitian matrix representing 
the effective Hamiltonian. In the following we shall choose units such that h = 1 and 
drop h. The requirement of completely positive evolution necessitates that the super- 
operators D(F m ,F n )p(t) be of the form [2011111 122 




(1) 



n , rn — 1 



®(F m ,F n )p(t) 



F n pFl-\{pFlF n + FlF n p) 



(2) 
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where {F n } n _[ is an orthonormal basis for all trace-zero operators on W, and the 
coefficient matrix (f mn ) is positive. The latter requirement implies f^ n — f nm , i.e., 
we can write f mn = a mn + ib mn , where (a mn ) is a symmetric real matrix and (b mn ) 
is a real-anti-symmetric matrix. 

Taking {<?k}k=i to be a basis for the Hermitian matrices on the Hilbert space % 
the master equation (JlJ with dissipation term |2| can be written in coordinate form 
as a linear matrix differential equation (DE) 23 

|f=(L + D)r, (3) 

where r = (r„) £ R n2 with r n = Tr(<r n p) and L and D are N 2 x N 2 (real) matrices 
with entries 

L mn = Tr(iH[a m ,a n }) = ^h n > Tr(ia n >[a m , a n ]), (4a) 

n' 

D mn = ^ fm'n' Tr(F^, cr m F n /CT„ - \ F^, F n > {a m , a n }) . (4b) 

Here {^4, B} = AB + BA and [A, B] = AB — BA are the usual matrix anti-commutator 
and commutator, respectively. If we choose a basis such that cr^ra = an< i the 

remaining basis elements form a basis for the trace-zero Hermitian matrices then 

r N2 = Tr(p) = ^ (5) 

is constant, i.e., f^i = 0, and we can define a reduced so-called Bloch vector 
r = (ri, . . . , r^v2_ 1 ) T [24], and rewrite the linear matrix DE for r as affine-linear 
DE 

r(t) = Ar(t) + c, (6) 

where A is an (iV 2 — 1) x (./V 2 — 1) real matrix with A mn — L mn + D mn and 
c rn = -^D mN 2. We can without loss of generality choose F m = a n . 

The objective of complete identification of the dynamics is thus reduced to 
identifying the coefficients h n > and /,„'„' of the N 2 xN 2 (Hermitian) coefficient matrix. 
Positivity constraints restrict the parameter space further to a convex subset of M. N . 



3. Identifiability given Complete State Information 

The Bloch equation for dissipative quantum systems always has at least one steady 
state r ss |25j for which we have 

= r ss = Ar ss + c. (7) 
This implies c = — Ar ss and shows that the Bloch vector r(i) is given by 

r(t) = r ss + e* A (r(0)-r ss ). (8) 

We note that if there are two steady states ril' 2 ^ then AAr ss = A(r S s' ) — riV) = 
shows that Arss'' = Aris" 1 and thus c is well-defined and 

r(t) = ri 2 J+e tA (r(0)-ri?) 

= r« + Ar ss + e tA (r(0)-r«- Ar ss ) 
= r«+e M (r(0)-rW) 
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noting that Ax. — implies cxp(tA)x = lx. 

Generically the superoperator A has d = N 2 — 1 distinct eigenvalues X n with 
corresponding eigenvectors v„. The eigenvalues and eigenvectors may be complex 
but as A itself is real, complex eigenvalues must occur in complex conjugate pairs. 
Suppose there are p < d/2 pairs of complex eigenvalues A* = -f n ± iuj n and d—2p real 
eigenvalues \ n = 7 n for n = p + 1, . . . , d — p. Writing the corresponding eigenvectors 
as v± = x„ ± iy n for n G h = {1, . . . ,p} and v n = x„ for n G 1% = {p + 1, • • • , d-p} 
with x„, y„ G M. d , and noting that A* — A shows that 

A(x + iy) = (7 + ico)(x + iy) =>■ A(x — iy) = (7 — iw)(x — iy). 

The vectors {x„}^~^ and {y n } p n=1 together form a (non-orthogonal) complete basis 
for R d and we can expand the initial state r and steady state r ss with respect to it 

d—p p d—p p 

r 8s = ^2p n x n + ^2/3' n y n , r = ^(a n + /3„)x„ + ^(a'„ + P' n )y n , (9) 

n— 1 n— 1 n— 1 n— 1 

where a„ and /3 n are real coefficients. This yields 



= r - r ss = ^2 a « x » + + E a « x « 

riG/i 116/2 



V" ^r v + + v -i + ^r v +_ v - 

/ t 2 L n 1 nJ 1 2i 
riG/i n£/ 2 



riG/i 116/2 



E 



where = \{a n =F ia^)v^ and v n = a„v„ and = if a„ — a' n = and £ n = 1 
otherwise. Inserting the last expression into ([8j and recalling Av^ = ("fn ± «w„)v^ 
we obtain 

p d—p 

r(i) =a + £ n e 7 ™* [a„ cos(u;„i) + b„ sin(w„t)] + £„a„e 7ret (10) 

n— 1 n— p+1 



where the coefficient vectors are ao = r ss , a„ = 2x„, b„ = — 2y n for n G ii, a„ = v 
for n G I 2 , and v± = x„ ± iy„ with x n = §(a n x n + a4y„), y„ = i(-a^x„ + a„y„) 
The coefficient vectors ao, a„ and b„ can be estimated along with the parameters u> 
and 7„ using Bayesian estimation, choosing the basis functions 



n 



n 



9o = 1 






(Ha) 


92n-i(t) = e 7 "* sm(u) n t) 


n = I, . 


■ • ,P 


(lib) 


92n(t) = e 7 "*cos(w„i) 


n = 1, . 


■ ■ ,P 


(11c) 


92p+n(t) = e 7 "*, 


n = 1, . 


..,N-1. 


(lid) 


Thus, provided 7^ for all n, we 


can determine all of the eigenvectors 


v„ and 


eigenvalues of A, so that A = SDS -1 


with 5 = , v 1 , . . . 


■v+,v p ,V p+ i,.. 


• > Vp+d) 
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being a matrix whose columns are the complex eigenvalues of A and D = 
diag(A+, X[ 



. A+, A~, 7 P +i, • • ■ , 7p+d) is a diagonal matrix with the corresponding 
eigenvalues A* = j n ± iui n . S will be invertible provided £„ 7^ for all n so that we 
have determined all eigenvectors. 

Theorem 1 The trajectory r(t) of an initial state uniquely determines the A and 
c = — Aslq provided r(0) — r ss , where r ss is a steady state, has non-zero overlap with 
all eigenvectors of A. 

If A has fewer than TV 2 — 1 distinct eigenvalues then A(t) = r(t) — r ss still 
satisfies a homogeneous linear equation A(t) = AA.(t) and we have A = SJS^ 1 , 
where S is a matrix determined by the (generalised) eigenvectors but J — diag( J n ) is 
the Jordan normal form of A consisting of irreducible Jordan blocks J„ of dimension 
k n with some eigenvalue A„. In this case we have A(i) = Se Jt S -1 A.(0), where S is a 
matrix whose columns are generalised eigenvectors v n k satisfying (A — X^I) k v nk = 0, 
(A — I) k ~ lv n k 7^ 0. The generalised eigenvectors can be chosen such as to be span 
l w _1 . The matrix exponential e Jt is block-diagonal with fc„ dimensional blocks 



£„(*) 



i?„ 




L J.2 D 
2 ^ -'"-r) 



4 2 J 
tR n 

Rn 



1+3 



tR n 



(12) 



where 



i?,„ = 1 for uj n = 0, 



cos(w„i) — sin(w n t) 
sin(w rl i) cos(w n t) 



for u> n =/= 0. 



(13) 



For non-generic A there are many different subcases that arise even in the simplest 
case of a single qubit. These are discussed in detail in Appendix A Most of these 



special cases rarely arise in practice and the nature of the positivity constraints impose 
further restrictions. Physically allowed examples of non-diagonalisablc Bloch matrices 
do exist. For qubits they involve extremely high decoherence rates leading to over- 
damped dynamics 23.1 , which are usually not of interest for quantum information. In 
higher dimensions there are more possibilities for A to non-diagonalisable although 
generically A is still diagonalisable. The analysis in the appendix suggests that the 
trajectory r(t) of one initial state still uniquely determines A and c = — Aslq even in 
the non-generic case provided (i) A does not have two or more generalised eigenvectors 
of the same degree k corresponding to the same eigenvalue A, and (ii) r(0) — r ss , where 
r ss is a steady state, has non-zero overlap with all (generalised) eigenvectors of A. 



4. Identifiability given Partial State Information 

A more interesting question is how much information about the dynamical generators 
we can obtain from partial state information such as measurement of a single 
observable. To address this problem it is useful to consider the identification problem 
in the Laplace domain |10) . The Laplace Transform of the Bloch equation ^ gives 

s 2 R(s) - sr(0) = AR(s)s + c = T ( sK j s A (14) 
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where R(s) is the Laplace transform of r(t) and T — (A, c) is the matrix A with the 
column vector c appended. Thus, given R(s) and r(0) then T and thus A and c are 
in principle fully determined. Indeed we only need to know R(s) for d + 1 values of 
Sk in general. Define the d + 1 x d + 1 matrices 

c , = /siR(si) ... s d+ iR(s d+ i)A ^ 15 ^ 

D = (sfR(si) - sir(O), ..., S 2 +1 R( Sd+1 )- Sd+ir (0)). (15b) 

We have D = TC and if det(C) ^ then T = DC^ 1 and thus we have fully determined 
both A and c. Thus if we could measure in the Laplace domain rather than the time 
domain only d + 1 measurements would be required to completely determine the 
dynamics although it not easy to see how to perform such measurements in practice. 



But ( 14 ) can also be useful if we have incomplete information about the state, 



e.g., if we can only measure some components of the Bloch vector. We can use ( 14 1 to 
express the unknown components in terms of the known ones. For example, suppose 
we know only the last component R d (s). Let R'(s), r'(0) and c' denote the first d — 1 
components of R(s), r(0) and c, and partition the matrix A as follows 



Then ( 14 ) can be rewritten as 

sR'(s) - r'(0) = A'TL'(s) + aR d (s) + c'/s (17a) 
sR d {s) - r d (0) = b T R'(s) + a dd R d {s) + c d /s. (17b) 

If s is not an eigenvalue of A' then the first equation gives 

R'( s ) = ( s /-^')" 1 [r'(0) + aii d (s) + c7s] (18) 

and inserting this into the second equation yields a non-linear equation in s that 
depends only on R d (s), 

= b T (sI - V(0) + ssLR d (s) + c'} - s 2 R d (s) + sr d (0) + sa dd R d (s) + c d , (19) 

which can be rewritten as 

p , x b^s/-ArV'(0) + c']+^(0)+c d 

Rd{s) ' S 2_ a(b T (sJ _ j40 -i a + add) • ( 2 °) 

The equation must be satisfied for all s. R d (s) is a rational function and we can 
therefore expand the numerator and denominator in terms of powers of s and match 
the coefficients, leading to a set of algebraic equations for the coefficients. We can 



identify which parameters of A contribute to the RHS of Eq.(20), hence determine 
what information can be extracted in principle. 

In practice we can use knowledge about A to infer the structure of the expected 
signal. From Sec. [3] we know that in the generic case A has iV 2 — 1 distinct eigenvalues. 
More precisely, we expect N— 1 real eigenvalues and (N 2 —N—2) /2 = N(N— 1) /2— 1 =: 
p pairs of complex conjugate eigenvalues and thus a signal r(t) of the form 



N-l 



r d (t) = c + ^2 c„e 7 "* + ^ e 7 "*[a n cos(uj n t) + b n sin(w„i)]. (21) 



n—1 n—1 
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Given stroboscopic measurement data rd(t k ) we can now use parameter estimation 
techniques to estimate oj n , j n and j' n and the coefficients a n , b n and c n . 

The approach above can be generalised to arbitrary partitioning of the observed 
signals and in principle we could consider more general observables. 



5. Application: Identification of Qubit Dynamics 

We illustrate the ideas above by applying them to a qubit evolving under general 
Markovian dynamics. Taking 

;;)•».-(; 7) ■ »-(; -°o <-> 

to be the unnormalised Pauli matrices, we can expand the Hamiltonian H = h^li + 
h-x&x + hyljy + h z a z , where hj(j — 0,x,y,z) are real. Choosing F n to be dissipation 
operators to be the normalised Pauli matrices, F\ = ~^o~x, F% — ^o~y> F*3 = ~^°~z> 
F 4 = -^h, and setting f mn = f£ n + if^ m , where (/^„) is a real symmetric matrix 
and (fmn) are real-anti-symmetric matrix, the general Bloch operators take the form 

A = /* + 2h z -(/* + /*) /* - 2^ , c = n/2 -/f 3 . (23) 
V /il - 2/1, /i + 2^ -(/* + /*)/ V // 2 / 

The anti-symmetric part of ^4 determines the Hamiltonian while the symmetric part of 
A determines the real parts of the dissipation coefficients f mn and the inhomogeneous 
part determines the imaginary part of f mn . 

Generically, we expect one real eigenvalue — 6 and a pair of complex eigenvalues 
—7 ± icu with <5, 7 > 0. Therefore we choose the basis functions for the Bayesian 
estimation 

S0(t) = l, gi(t) = e-i* cos(wt), g 2 (t) = e-T t sin(wt), g 3 (i) = e" 5t , (24) 

i.e., in particular the signal should be a linear combination of these basis functions 
s(t) = aog (t) + ai<7i(t) + CL292(t) + a 393(t)i an d the Laplace transform of the signal 
should be the corresponding linear combination of 

G (s) = Gi(s) = - — S ^ 2 7 — 2 , G 2 (s) = ^ — G 3 (s) = — j-r. 

The parameter estimation in all of the following cases uses a standard Bayesian 
estimation technique which involves orthogonal projections of the data onto the basis 
functions and optimisation to maximise the log-likelihood, this method being described 
more fully elsewhere (TTJ [12] . 



5.1. Identification given full Bloch vector information 

We first consider a generic qubit system for which the Bloch matrix A has three 
distinct eigenvalues. The system is prepared in a fixed initial state ro and allowed 
to evolve under the unknown dynamics. After evolving for different times, t k , the 
components of the Bloch vector ri(i fe ) = (a x (t k )), r 2 (t k ) = (a y (t k )), r 3 (t k ) = (a z (t k )) 
are measured. For each time point the experiment is repeated N e times to determine 
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t (arb. units) 

Figure 1: Model 1, T = 50, N = 1000 and N e = 1000. Single input state but full 
tomographic data. For longer sampling times, the signal to noise ratio is too low and 
reconstruction accuracy suffers. 



the relative frequency of measurement outcomes and 1. If we can measure all three 
components of the Bloch vector (equivalent to QST of a qubit) then we obtain three 
data traces of a form shown in Figure [T] 

From these data traces, we then estimate the underlying signal parameters and 



express the signal in the form ( 10 1 and consequently reconstruct the matrix A and c 



using the technique outlined in Section [3j Figure [2] shows the cumulative probability 
density of the relative root-mean-squared error of the reconstructed Bloch matrix A 
and constant c for a typical generic system, specifically (Model 1) 



'-0.0650 -2.0000 2.0300 \ 
1=1 2.0000 -0.0650 -4.0000 , 
^1.9700 4.0000 -0.0900/ 



c = 



-0.0424N 


0.0636 




for different choices of the signal length, number of samples N and experiment 
repetitions N e , A set of 100 randomly generated data traces were used for 
the reconstruction. Time samples for each run are chosen by low-discrepancy 
sampling [26 . As the figure shows the median error for N = N e = 1000 is 1.5%; 
for N = 1000, N e = 5000 the relative error of more than 95% of all runs is below 
1%. There is an optimal range of the signal length around T = 50 for this system. 
The reconstruction error increases for substantially shorter or longer signals; in the 
latter case this is due to the signal vanishing. Comparing the results for T = 50 and 
N = 2000, N e = 1000 and N = 1000, N e = 2000 shows that there is virtually no 
difference in the probability density, suggesting that it makes little difference if we 
increase the number of time samples or the signal to noise ratio of the sample points. 



Open quantum system identification 



10 



100 



80 



c 

a; 



-Q 

-Q 

O 



60 



40 



20 



L 
10" 




N500, Ne500 
N1000, Ne1000 
N1000, Ne2000 
N1000, Ne5000 
N1000, Ne50000 
N250, Ne1000 
N500, Ne1000 
N500, Ne500 
N500, Ne500 



10' 



10 10 
relative error (%) 

Figure 2: Reconstruction-error cumulative probability density for model 1. We 
examined different sampling times, sampling points and experiment repetitions. As 
expected, as the noise of each data point decreases, and the number of time samples 
increases, the reconstruction accuracy diminishes. Different system models give similar 
error probability plots. 



5.2. Identification based on single component of Block vector 

Suppose we can only measure one component of the Bloch vector in a fixed 
measurement basis. This is typically the case where there is a physically preferred or 
engineered mechanism and before coherent rotations of the state are possible. Without 
loss of generality assume we measure rd(t) = r z (t), for instance by atomic population 
measurement by fluorescence shelving [27] . We can write the Laplace transform of the 
signal as 

r M — — + fl i( s + 7) + , fl 3 _ C 3 s 3 + C 2 s 2 + Cis + Cp 
K d{s)- s + (s + 7)2+w 2 + s + 6 - a* + D 3 s 3 + D 2 s 2 + D l3 ' ( ' 

Assuming we have estimated the coefficients clq, a\, 03 as well as u>, 7 and S using 
signal parameter estimation, the constants Ck, -Dfc for k = 0, 1, 2 are determined from 
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the observed signal 

C = a S(j 2 + uj 2 ) (26a) 

C\ = a ( 7 2 + 2j5 + ui 2 ) + ai 5j + a 2 Suj + a 3 ( 7 2 + uj 2 ) (26b) 

C 2 = a (S + 27) + 01(7 + 5) + a 2 uj + a 3 2-f (26c) 

C 3 = a + ax + a 3 (26d) 

D 1 = 5{ 1 2 +uj 2 ) (26e) 

L> 2 = 2 7 <S + 7 2 + cj 2 (26f) 

^3 = 2 7 + 5 (26g) 



Comparing with (201 for a generic single-qubit Bloch matrix A — (a mn ) shows that 



the constants must satisfy 

C =(a 2 ia 32 - a 22 a 3 i)ci + (a 12 a 3 i - a n a 32 )c 2 + (a n a 22 - ai 2 a 21 )c 3 (27a) 
Ci =a 3 ici + a 32 c 2 - (an + a 22 )c 3 + (a 32 a 2i - a 3 ia 22 )r x (0) 

+ (031012 - a 32 an)r y (Q) + (ana 22 - 012021)^(0) (27b) 

C 2 =c 3 + a 3 ir K (0) + a 32 r y (Q) - (an + a 22 )r z (0) (27c) 
Di =a 13 a 3 ia 22 - a i2 a 23 a 3 i - ai 3 a 2i a 32 + ana 23 a 32 + ai 2 a 2i a 33 - ana 22 a 33 ( 27d) 

D 2 = - a i2 a 2 i - a i3 a 3 i - a 23 a 32 + a 22 a 33 + an (a 22 + a 33 ) (27c) 

D 3 = - an - a 22 - a 33 (27f) 

As the constants C&, Dfc and the initial state r x (0), r y (0) and r^O) are known, we have 
a system of six non-linear algebraic equations for the coefficients of A and c. Given 
that A in general has 9 coefficients and c has three, these equations do not determine A 
uniquely, although we may be able to effectively identify all unknown parameters if we 
have prior information about the dominant types of dissipation, effectively restricting 
the form of A and c. This may be derived from physically motivated grounds or from 
prior measurements of similar systems. We apply this to some examples below. 



5.2.1. Dephasing in a general basis Consider the special case of a qubit subject to a 
fixed but unknown Hamiltonian H = h{h x a x + h y a y + h z a z ) and dephasing operator 
V = -^{&a x + fl&y + l&z), which corresponds to dephasing in an arbitrary basis. In 
this case the master equation can be written in Lindblad form 

£p(t) =-{[H,p] + Vp(t)V^ - \{V^V, P (t)} (28) 

and expanding with respect to the normalised Pauli matrices, the corresponding terms 
in the Bloch equation are 

/-(/3 2 +7 2 ) -h z +a/3 h y + a-f \ 
A=\ h z +af3 -(q 2 + 7 2 ) -h x + £7 , c = 0. (29) 
\-h y + aj h x +fij -(a 2 + f3 2 )J 



i.e., we have a total of six parameters. We see from (271 that for c = there are only 
five non-trivial signal parameters. Thus, a single trace is not sufficient to completely 
determine all of the model parameters. If the system is initialised in a measurement 
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Figure 3: Comparison of cumulative probability density for the relative Hilbert- 
Schmidt error of estimated Bloch matrix A for Model 2. We examined different 
starting states and reconstructions based on full or partial observable information. 
The reconstruction based on full XYZ information did not assume any particular 
decoherence model, unlike the reconstruction based on only partial observables. Initial 
states were either ro = (1, 0, 0) T or ro = (0, 0, 1) T . Partial observation was either the 
z trace only, y and z, x and z, or x and y traces. 



basis state r = (0,0, 1) T then (26) and (27) give explicitly 



Co=0 

Ci =fr 2 + 7 2 (a 2 + ^ 2 + 7 2 ) 
C 2 =a 2 + (3 2 + 2 7 2 

D\ =h 2 z (a 2 + f3 2 ) + h 2 y (a 2 + 7 2 ) + h 2 x {p 2 + 7 2 ) 

— 2h x h z a r y — 2h x h y af3 — 2h z h y /3 r y 
D 2 =h 2 x + h 2 y + h 2 z + (a 2 + (3 2 ) 2 + 7 2 (2a 2 + 2/3 2 + 7 2 ) 
D 3 =2(a 2 +/3 2 + 7 2 ) = 2(C 2 - 7 2 ) 



(30a) 
(30b) 
(30c) 

(30d) 
(30e) 
(30f) 



We cannot determine all 6 model parameters, the observed signal is invariant under a 
rotation around the z-axis, hence we cannot determine the azimuthal components of 
the Hamiltonian and dephasing simultaneously, only their relative orientation. This is 
similar to the original situation of Hamiltonian identification where there is a "gauge 
symmetry" which exists due to only being able to initialise and measure in a fixed 
basis, i.e. we do not have a phase reference to measure the x and y components of the 
dynamics separately. However, we can determine h 2 + h 2 and a 2 +(3 2 . We can "gauge 
fix" , e.g. h y = 0, and define a and j3 in reference to the projection of the Hamiltonian 
onto the x — y plane. 
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Figure 4: Comparison of probability density for the relative Hilbert-Schmidt norm 
error of estimated Bloch matrix A for Model 3. Reconstruction procedures as in 
Figure [3j The optimum sampling time is shorter than for Model 2 as the signal decays 
faster. 



We have solved for the analytic solutions of (30) but due to signal noise, these 
may not be exactly consistent. We therefore implemented the non-linear system of 
equations by minimising the least-squares errors using numerical optimisation starting 
with an initial guess derived from the analytic solutions. We plot the cumulative 
probability distribution for the errors in Figure [3] for Model 2, 



.4 



which corresponds to h x = h z = 2, h y — and a — 0.1, j3 — j — 0.2. We implemented 
various initial states, observables and prior knowledge. 




5.2.2. Dephasing and Relaxation in Measurement basis We also consider the special 
case of a qubit sub j ect to a fixed but unknown Hamiltonian H = ^(h x a x +h y a- y +h z a z ) , 
and relaxation and decoherence along the measurement basis governed by the Lindblad 
equation 

jp{t) = - l -[H,p] + Y i V k p{t)Vl \{Vlv k , P (t)} (31) 

k=l 

with dephasing operator V\ — VTfj z , and relaxation operators V2 = and 
V3 = \/72|l)(0| where we include both spontaneous emission and absorption. Again, 
expanding with respect to the normalised Pauli basis the Bloch operators take the 
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Figure 5: Error histogram plot (Model 3, r = (0,0, 1) T , z-trace, N — N e — 1000, 
T = 15). The absolute errors in the elements of matrix A were comparable with 
each other, but as the Hamiltonian parameters are much larger than the decoherence 
parameters, the relative errors of the latter are commensurately greater. Due to 
rotational symmetry around z, we can only determine \ h 2 . + h 2 . 



forr 



— r c ff 


-h, 


hy \ 




K 


— r c fr 


-h x 


c= f 


-ky 




-Is) 


\A 7 



(32) 



where T eB = -2T - ±A 7 , A 7 = - 72 ), 7s = 7l + 72 . For r = (0,0, 1) T the 

equations are 

Co =(r c 2 ff + h 2 z )A 1 (33a) 

d =2r cff A 7 + r 2 cS + h 2 z (33b) 

C 2 =A 7 + 2r cff (33c) 

D t = ls {hl + r c 2 ff ) + {hi + h 2 y )r eS (33d) 

D 2 =T 2 cS - 2r cff7s + hi + h 2 y + h 2 z (33e) 
D 3 =2r eff + 7s . (33f) 

Analytical solutions to these have also been found but as in the previous case, 
numerical optimisation is required due to signal noise. Again, we can only determine 
hi + h 2 and not h x and h y individually if we start off with Tq — (0, 0, 1) T and measure 



r z . If we start off with rrj = (1, 0, 0) , we obtain the same expressions as in (33 1 except 
for the two equations 

C* 2 =A 7 - h y (34a) 
Ci =2r cff A 7 + h x h z - h y T cS , (34b) 
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which allow us in principle to identify all six parameters. However, the equations 
appear to be sensitive to noise and often lead to incorrect model parameters. This 
behaviour is a subject for further research. The results for Model 3 

/-0.25 -3.00 2.0 \ / \ 

A = 3.00 -0.25 -1.0 , c = , 
\-2.00 -1.00 0.1 / \0.1/y/2/ 

are shown in Figure [4] for various initial states, observables, and prior knowledge. 
Figure [5] shows the histogram of the error distribution of the model parameters for 
one reconstruction setting. 



5.3. Identification given two components of the Block vector 

In some systems, e.g. bulk nuclear magnetic resonance, determination of the x and y 
components of the Bloch vector is easy to measure. In this case, we assume we can 
measure both r x (t) and r y (t). Inserting (18b) into (18a) gives 

sK'(s) - r'(0) = A'K'(s) + a[b T R'(s) + a dd R d (s) + c d /s + r d (0)} + c'/s. 

Solving for R'(s) gives 

R'(«) = [si' - A' - ab T /(s - a^rVW + c'/s + a(c d /s + r d (0))/(s - a dd ) (35) 

from which we can obtain two explicit rational equations for R x (s) and R y (s). The 
denominators of both are the same and the coefficients of the polynomials on the top 
and bottom can be related to the observed signals as before. 



For the two examples in sections 5.2.1 and 5.2.2 we implemented the 



reconstruction procedure, the results of which are presented in Figures [3] and [4] For 
the sake of comparison we considered not only x and y but different combination of two 
traces including y and z, and x and z observable information for the reconstruction. 
The general observation is that for the systems considered the z-trace appeared to 
be the best choice for the reconstruction. The addition of the extra information 
in the form of the x or y trace did not improve the reconstruction accuracy over 
reconstruction using the z trace alone for these examples, and using the x and y 
traces instead appear to result in substantially worse reconstruction results. There are 
also initial state effects with initialisation in the z-eigenstate being decided preferable 
compared to initialisation in the s-eigenstate in the second case, for example. 



6. Summary and Discussion 

We have studied the problem of identification of the dynamics of open systems in 
a Markovian environment without recourse to quantum process tomography. By 
utilising stroboscopic measurements over time we showed that the system dynamics 
is completely determined by the evolution of a single initial state, at least in the 
generic case when A has distinct eigenvalues and the initial state has overlap with 
all eigenspaces of A, i.e., we can in principle fully determine the Bloch operator A 
and inhomogeneous term c from stroboscopic measurements of the state at different 
times. Analysis of non-generic cases suggests that this result basically still holds in 
most cases and we can give specific conditions under which full information about 
the dynamics can be obtained. For systems governed by Markovian evolution this 
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approach is preferable to process tomography for several reasons. First, it does not 
require initialisation in a large set of different basis states. Second, it avoids the rather 
messy procedure of extracting the dynamical generators from process tomography 
data, which is susceptible noise. Third, it is more efficient than performing process 
tomography at numerous different times. 

By going to the Laplacian domain, we can identify algebraic equations relating the 
observed signals to the models parameters. In the absence of a mechanism of directly 
sampling in the Laplacian domain, we can use Bayesian signal estimation in the time 
domain to find a functional representation of the measurements more amenable to 
Laplacian domain analysis. The most important advantage of this approach is that it 
shows that in many cases most or all model parameters can be obtained from far less 
data, i.e., we do not even require full state tomography. Rather we can extract the 
same information from measurements of a smaller set of observables. This approach 
is especially useful if there is physical information restricting the type of decoherence 
model we expect, and the dynamical generators depend on fewer parameters. Our 
analysis of the qubit cases shows that not only is it possible to extract all of the system 
parameters from this restricted data but often this even leads to better reconstruction 
results than using full state reconstruction. 

In part these results can be explained by overfitting, i.e., trying to fit the most 
general model given noisy data when the actual dynamics is determined by a smaller 
set of parameters leads to fitting the noise and results in worse model over all. 
However, overfitting does not appear to explain the peculiar dependence of the results 
on the choice of initial state or observables noted earlier. For many of the systems 
considered the success of the reconstruction showed a strong dependence on the choice 
of the initial state and the type of measurement, with initialisation in a z-eigenstate 
of measurement of (o~ z ) being distinctly preferable over other combinations such as 
measuring (a x ) or initialising in an x-eigenstate. This effect was most pronounced 
for system subject to relaxation in the measurement basis. In this case the choice of 
the initial state (1,0, 0) T in principle breaks the rotational symmetry, allowing the 
determination of all six model parameters, instead of the five using the initial state 
(0,0, 1) T . However, in the former case the reconstruction often fails, converging to 
a model which is distinctly different from the actual model. There appears to be 
a trade-off between the breadth of the information one versus the quality of the 
information gathered but the issue of the optimal choice of the initial state and 
observables warrants further study, as do the choice of optimal signal lengths, adaptive 
sampling, stability of reconstruction algorithms in the presence of noise among others. 
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Appendix A. Non-diagonalisable Bloch matrices for N = 2 

As the complex eigenvalues of the Bloch matrix A must occur in pairs 7 ± iu>, we can 
only have a non-trivial Jordan block for a qubit system if uj = 0, i.e., the eigenvalues 
are real, and there are five cases: A = SJS^ 1 with 

'71 0\ h 1 0\ / 7 1 N 

Ji = I 7l 1, J 2 = I 7l , J 3 = 7 1 , J 4 = 7 




72/ \0 7/ \0 7; 

(A.l) 

and J5 = 7/. The corresponding time evolution operators are e tA — Se tJ S 1 with 
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; ; = I r'" 7 I . = I 

Case 1. There are two proper eigenvectors xn, x 2 and one generalized 
eigenvector x 12 satisfying (A — 7„/)x„ 1 = for n — 1,2 and (A — 7 1 /) 2 x 12 = 0. 
Expanding with respect to these A(0) = anxn +ai 2 xi 2 + a 2 x 2 and using the Jordan 
form of A above gives 

r(t) = r ss + e* 71 (a n + a 12 t)x n + e* 7l ai 2 x 12 + e* 72 a 2 x 2 
= r ss + e* 7l toi 2 xn + e t7l (aiixn + ai 2 xi 2 ) + e t72 a 2 x 2 
= a + te t7l ai + e t7l a 2 + e* 72 a 3 



Using parameter estimation we can in principle estimate 7„, n — 1,2, as well as the 
coefficient vectors a„ for n = 0, 1, 2, 3. a determines the steady state and ai and a 3 
determine the two primary eigenvectors as we can absorb the factors a i2 and a 2 - a 2 
is a generalised eigenvector as (A — 7i/)(aiiXii + ai 2 xi 2 ) = (A — 7i/)ai 2 xi 2 7^ 0, 
(A — 7i/) 2 (aiixii + ai 2 xi 2 ) = 0. Thus we have S = (ai, a 2 , a 3 ) and both eigenvalues 
and A is completely determined. 

Case 2. There are two proper eigenvectors xi, x^ with eigenvalue 71 and one 
proper eigenvector x 2 with eigenvalue 7 2 . Expanding with respect to these eigenvectors 
A(0) = ctixi + o^Xi + a 2 x 2 and 

r(t) = r ss + e* 71 (aiXi + a'^) + e t72 a 2 x 2 
= a + e* 7l ai +e t72 a 2 

Using parameter estimation we can estimate 7„ as well as the coefficient vectors a„ 
for n — 0, 1,2. a determines one steady state and a! and a 2 determine two of the 
proper eigenvectors of A but without further information we cannot determine the 
third eigenvector. 

Case 3. In this case there is only one eigenvalue 7 and one proper eigenvector xn 
and two generalised eigenvectors x i2 , xi 3 . Expanding with respect to these generalised 
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eigenvectors, A(0) = anXn + ai 2 xi 2 + ai 3 Xi 3 , and using the Jordan form gives 

r(t) = r ss + e* 7 (a 11 + a 12 t + \a ld ,t 2 )^n + e* 7 (a 12 + tai 3 )x 12 + e* 7 a 13 x 13 

= r ss + iai 3 t 2 e t7 xn + ie* 7 (ai 2 xn + ai 3 xi 2 ) + e' 7 (aiixii + a 12 x 12 + ai 3 xi 3 ) 
= a + tV 7 ai + te t7 a 2 + e t7 a 3 

Using parameter estimation we can estimate 7 as well as the coefficient vectors a„ 
for n = 0,1,2,3. Again ao determines a steady state, ai = ^o;i 3 Xii is a proper 
eigenvector of A, a 2 and ai are generalised eigenvectors with (A — -fl) 2 a 2 = and 
(A — 7/) 3 a 3 = 0, respectively. S — (ai, a 2 , a 3 ) and 7 and A is completely determined. 

Case 4. In this case there is again only one eigenvalue 7 but there are two 
proper eigenvectors xi, and one generalised eigenvector x i2 with (A — ■yl)xi 2 7^ 
but (A — -fl) 2 xi 2 = 0. Expanding with respect to these generalised eigenvectors, 
A(0) = ctixi + a'jx'j + ai 2 xi, and using the Jordan form gives 

r(t) = r ss + e* 7 [(ai + ai 2 i)xi + ai 2 xi 2 + a'^] 

= r ss + fe t7 ai 2 xi + e* 7 (aiXi + a[x[ + ai 2 xi 2 ) 
= a + te* 7 ai + e* 7 a 2 

Using parameter estimation we can estimate 7 as well as the coefficient vectors a„ for 
n = 0, 1, 2. Again a determines a steady state, a! = iai 3 x n is a proper eigenvector 
of A, a 2 determines another generalised eigenvector but the third eigenvector can 
again not be determined. 

Case 5. In this case any vector is an eigenvector with eigenvalue 7 and we have 
r(t) = r ss + e* 7 A(0). Using parameter estimation we can determine 7, a and A(0) 
and in this case this information is sufficient to determine A and c — if it is known 
that we have case 5. 

Cases 1 and 3 can clearly be identified from the structure of the observed signal. 
A signal of the form observed in case 2 can also arise if A has three distinct (real) 
eigenvalues but the initial state is such that A(0) has no overlap with one of the 
eigenvectors. Case 4 can also arise if A has two real eigenvalues, one of which 
corresponding to a Jordan block of size 2 but the initial state is such that A(0) has 
no overlap with one of the proper eigenvectors. Finally, a signal of the form observed 
in Case 5 can arise not only when A is a scalar matrix but also when A(0) docs not 
have overlap with all (generalised) eigenvectors of A. 
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